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We investigate the dynamics of photoexcited carriers and nonequilibrium phonons in graphene 
by solving the microscopic kinetic Bloch equations. The pump and drift effects from the laser field 
as well as the relevant scatterings (including the Coulomb scattering with dynamic screening) are 
explicitly included. When the pump-photon energy is high enough, the influence of the drift term 
is shown to be negligible and the isotropic hot-electron Fermi distribution is established under the 
scattering during the linearly polarized laser pulse investigated here. However, in the case with 
low pump-phonon energy, the drift term is important and leads to a net momentum transfer from 
the electric field to electrons. Due to this net momentum and the dominant Coulomb scattering, a 
drifted Fermi distribution different from the one established under static electric field is found to 
be established in several hundred femtoseconds. We also show that the Auger process investigated 
in the literature involving only the diagonal terms of density matrices is forbidden by the dynamic 
screening. However, we propose an Auger process involving the interband coherence and show that 
it contributes to the dynamics of carriers when the pump-photon energy is low. In addition, the 
anisotropically momentum-resolved hot-phonon temperatures due to the linearly polarized light are 
also investigated, with the underlying physics revealed. 

PACS numbers: 78.67. Wj, 78.47.J-, 71.10.-w, 72.80.Vp 



I. INTRODUCTION 

Graphene is a promising material which has received 
intensive attention both theoretically and experimentally 
due to its excellent transport and optical properties 
Especially, its linear band structure and zero band gap 
lead to unique optical and electrical properties, making 
it suitable for various optoelectronic applications >i£ A 
comprehensive understanding of the optical and electri- 
cal properties is one of the prerequisites for its potential 
applications. 

Experimentally, the time-resolved pump-probe mea- 
surement is widely used to investigate the dynamics 
of carriers in grapheneJ^r— With this method, both 
the differential reflection and the differential transmis- 
sion (DT) after pump are measuredJ^r— In the ex- 
periments with medium photoexcited electron density 
(around 10 12 cm 2 ) and probe-photon energy much 
higher than the equilibrium Fermi energy, the posi- 
tive DT is often observed with its fast relaxation of 
several hundred femtoseconds followed by a slower pi- 
cosecond relaxationj 13 i 20 i 21 i 25 The mechanism leading to 
these two distinct relaxations is attributed to the rapid 
carrier-phonon thermalization and the slow decay of hot 
phonons, respectivel y 20 ' 22 ^ 1 On the other hand, in re- 
cent years, experiments with higher photoexcited elec- 
tron densities (JVf > 10 13 C m- 2 )2iH show that a 
negative DT appears when the probe-photon energy 
(> 1500 meV) is much higher than the equilibrium 
Fermi energ y 23 ' 26 Such negative DT was usually claimed 
to arise from the renormalization of the single particle 
energy by the electron-electron Coulomb Hartree-Fock 



contributionj 23 ' 26 It is noted that although the observed 
DT with high probe-photon energy and low photoex- 
cited electron density is usually positive, negative DT 
was also reported with Nf! = 3 x 10 11 cm -2 by Plo- 
chocka et al. where the sample is 70-100 layer epitaxial 
graphene exhibiting the Dirac-type electron spectrum^ 
In contrast to the negative DT with probe-photon en- 
ergy much higher than equilibrium Fermi energy, other 
types of negative DTs with lower probe-photon energy 
have already been observed in previous works i i 15 Sun 
et al. reported the appearance of the negative DT un- 
der weak pump intensity with the probe-photon energy 
between 528 to 697 meV,— which is lower than twice of 
the Fermi energy (Ep ~ 350 meV). This negative DT 
was suggested to come from the weakening of the Pauli 
blocking due to the heating of the electrons by the pump 
pulsed However, this claim has not yet been verified by 
the microscopic theoretical approach. Soon after that, 
the negative DT was also observed by George et al. in the 
case when the probe-photon energy is as low as a few tens 
of meV (with the frequency about several terahertz). 15 
With such low probe-photon energy, the increase of the 
intraband absorption is assumed to be responsible for the 
negative DTJ£ In addition to these phenomena, with the 
pump-probe method, the optical gain has also been re- 
ported in optically pumped graphene very recently^ 

In addition to these experimental works, many the- 
oretical studies have been devoted to get a clear com- 
prehension on the dynamics of photoexcited carriers and 
phonons i 20 ' 25 ' 29 ' 31 " — Among these works, the simple rate 
equation approac h 20 ' 25 ' 39 assumes the establishment of 
the Fermi distribution as a starting point and, therefore, 
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fails to incorporate the pump process and the buildup of 
the Fermi distribution. The widely used time-dependent 
Boltzmann equation approac h 31 ! 34 ' 38 can fulfill the above 
aims but neglects the interband coherence^ which 
should be very important during the pump process. More 
elaborative microscopic investigations with the interband 
coherence included have been carried out based on many- 
particle density-matrix approach i 33 ' 35 " — However, these 
theoretical investigations are still far away from comple- 
tion, as stated in the following. 

Among those papers with the interband coherence 
included ) 33 ' 35 " — the scattering terms under the Marko- 
vian approximation with the Coulomb screening in the 
static limit are applied. Nevertheless, with the static 
screening, the electron-electron Coulomb scattering di- 
verges when the initial and final states of the scattered 
electrons satisfy the condition ?w_F|q| = |fax>| t 40 ' 41 where 
q and Huj are the differences in momentum and energy 
between the initial and final states, respectively. To avoid 
this divergence, in these papers, the delta function origi- 
nated from the Markovian approximation and reflecting 
the energy conservation in the Coulomb scattering term, 
was replaced by a Lorentzian function with the broaden- 
ing attributed to the non-Markovian effect^ 2 - Neverthe- 
less, under the more natural dynamic screening ) 31 ' 43 " — 
this broadening is not needed. This is because that the 
dynamic screening also diverges in a higher order when 
ftu_F|q| = |fio>| ) 46 ' 47 and therefore reduces the original di- 
vergent scattering amplitudes to zero. 

Another important issue frequently investigated is the 
influence of the Auger process. In previous works, the in- 
terband coherence was neglected and the Auger process 
was defined as the process that one of the two scattered 
electrons is transferred from one band to the other while 
the other electron remains in the same bandi 33 ' 34 ' 36 In 
the works by Winzer et aL ) 33 ' 36 such process was partic- 
ularly stressed. However, under the dynamic screening, 
this Auger process is totally forbidden because it also 
satisfies the strong restriction ^_p|q| = \hu>\ due to the 
energy and momentum conservation. However, in the 
case with the optical transition, the interband coherence 
is generated. With the interband coherence, there are 
Coulomb scattering processes with three out of the four 
band indices in the Coulomb interaction Hamiltonian be- 
ing the same and the left one being different. This is also 
a kind of Auger process. However, the influence of such 
process has not been addressed in the literature. 

Finally, in the microscopic equations obtained based on 
the many-particle density-matrix framework by Knorr et 
«L ) 33 ' 35 ~ — the pump term was obtained under the vector 
potential gauge. However, the drift term which describes 
the acceleration of the electrons by the electric field is not 
included. The drift term may be unimportant when the 
pump-photon energy is very high as investigated in these 
works. Nevertheless, it can be expected to be important 
when the pump-photon energy is low. On the other hand, 
in the time-dependent Boltzmann equation approach, the 
drift term is included but the pump term in the equa- 



tion has to be introduced based on a Fermi-golden-rule- 
like approximation and hence is suitable only for a weak 
pump intensity. 3 - Moreover, the Rabi flopping 4 ^ is absent 
in this approach. Therefore, a general equation based on 
the microscopic many-particle density-matrix approach 
that includes both the drift and pump terms naturally is 
needed. 

In this paper, we perform a microscopic investigation 
on the dynamics of photoexcited carriers and phonons 
in graphene by setting up and solving the kinetic Bloch 
equations with the electron-phonon, electron-impurity 
and electron-electron Coulomb scatterings explicitly 
included ! 48 ' 49 In our study, the dynamic screenin g 46 ' 50 " — 
is adopted in the Coulomb scattering. Moreover, with 
the gauge invariant approach,— both the drift and pump 
terms are obtained naturally. We look into the dy- 
namics of carriers photoexcited by the linear polarized 
laser pulses with the pump-photon energy in both near- 
infrared and THz regimes under medium pump intensity. 
When the pump-photon energy is high enough, the influ- 
ence of the drift term is shown to be negligible compared 
to that from the pump process. Moreover, the anisotropic 
photoexcited electrons tend to be isotropic under the 
scattering and an isotropic hot-electron Fermi distribu- 
tion is established before the end of the pulse investigated 
here. On the contrary, when the pump-photon energy is 
low, the drift term causes a net momentum transfer from 
the electric field to electrons. Together with the dom- 
inant Coulomb scattering, a drifted Fermi distribution 
is established in several hundred femtoseconds. More- 
over, the form of the drifted Fermi distribution differs 
from the one obtained under static electrics field^ 5 - Be- 
sides, we also show that the temporal evolution of the DT 
measured by Hale et al. can be well fitted with our mi- 
croscopic calculation. 25 In addition, although the Auger 
process involving only the diagonal terms of the density 
matrices is forbidden by the dynamic screening, we show 
that the Auger process involving the interband coherence 
still contributes to the dynamics of carriers. However, its 
contribution is important only when the pump-photon 
energy is low. Similarly, it is shown that the terms ne- 
glected in the rotation-wave approximation^ 8 - which is 
widely accepted in semiconductor optics, also become 
important in the case of low pump energy. We also in- 
vestigate the dynamics of phonons. Due to the linearly 
polarized pump pulse, the q-resolved hot-phonon tem- 
peratures are anisotropic. Moreover, the anisotropic hot- 
phonon temperatures under high pump energy are very 
different from those under low pump energy. Finally, we 
investigate the negative DT by fitting the temporal evo- 
lution of DT in the experimental work by Sun et al.^ 
Our results support their suggestion that the negative 
DT comes from the weakening of the Pauli blocking due 
to the heating of the electrons by the pump pulse. 

This paper is organized as follows. In Sec. II, we set up 
the model, lay out the kinetic Bloch equations and then 
give a simple analysis on the kinetic Bloch equations. In 
Sec. Ill the results obtained numerically from the kinetic 



3 



Bloch equations are presented. We summarize in Sec. IV. 



II. MODEL AND FORMALISM 

We start our investigation from graphene on Si02 and 
SiC substrates. In the absence of external field, the en- 
ergy dispersion of graphene near the K and K' points 
can be described by the effective Hamiltonian 53 (H = 1) 



ff^' AB (k) = v F (fiT x k x +T y ky) 



(1) 



Here, v F represents the Fermi velocity, [i = 1 (—1) 
for K (K r ) valley and r are the Pauli matrices in the 
pseudospin space formed by the A and B sublattices 
of the honeycomb lattice. The eigenvalues of H[f' AB 
are e^k = r)v F k and the corresponding cigenstates are 
Kr, = (w e "' A . l) T /\/2 with j] = 1 (-1) for conduction 
(valence) band and #k representing the polar angle of k. 
The additional electron-optical-filed interaction Hamil- 
tonian under the vector potential gauge is given by^ 



H, 



)i,AB 
photon 



(M) = \e\v F \p,T x A x {t)+TyA y (t)]/c. (2) 



Note that this Hamiltonian is consistent with the one 
obtained from the tight-binding approach j 55 ! 56 Here, we 
only investigate the linearly polarized light with the elec- 
tric field 



E(t) = -EoCos(w L t)exp[-t 2 /(2of)] 



(3) 



along the x direction. For simplicity, the vector potential 



is taken to be A(t) 



-Eosin(a; L t)exp[-t 2 /(2of)]. Our 



investigation is performed in the base set of the eigen- 
states of the Hamiltonian [Eq. (JTJ]. Under this base set, 
the Hamiltonians are given by 

Hg(\s.)=v F ka z , (4) 
#photon( k ) = \e\v F [A(t) ■ ka z + M k x A(t) • za y ]/c, (5) 

with <r representing the Pauli matrices in the conduction 
and valence band space, k indicating the direction of k 
and z standing for the unit vector along the z direction. 

Exploiting the nonequilibrium Green's function ap- 
proach with the gradient expansion as well as the gen- 
eralized Kadanoff-Baym ansatz, the gauge invariant ki- 
netic equations^ 8 - of the carriers under the Markovian ap- 
proximation rea d 45 ' 49 ' 57 ' 58 (see Appendix A for detailed 
derivation of the coherent and drift terms) 



coh 



drift 



(6) 



in which p M k stand for the density matrices with their 
diagonal terms p^k,!)!) = ffik.r/ representing the elec- 
tron distribution functions and the off-diagonal terms 
/0/^k.i,-i = P*^k _i i standing for the interband coherence. 
The drift terms can be written as 



dtPfj.k 



drift 



|e|E- V k p M k, 



(7) 



which describe the acceleration of electrons by the laser 
field E. The coherent terms are given by 



<9tPVk| con 
Here, [A, B] 



-i[v F ka z 



%k 



^Pump' P^k] 



(8) 



AB — BA is the commutator 
and flp ump = — \e\v F /jA x sin 9ua y /c comes from the 
laser field describing the pump process. = 

- E q Sk.k-qPMk-q^k-q.k^fa. °) stand for thc Hartree- 
Fock self-energies^ in which S^ k q are the form fac- 
tors with their elements being S^.k-q.^' = '^kJV'k-qr;' 
and V r (q, u>) = V^/e(q, w) represent the screened two- 
dimensional Coulomb potentials with V® = 2itv F r s /q be- 
ing the bare Coulomb potential. Here r s stands for the 
dimensionless Wigner-Seitz radiusi 51 ' 60 " — e(q, ui) is the 
dielectric function under the dynamic screening) 48 ' 51 i 60 
which is given b y 46 ' 50 " — 



e(q,u,) = l- £ Kk +q , 



2(//ik,r/ //jk+q,^')^? 



£>)k - ETj'k+q + LO + iO+ 

(9) 

It is noted that the laser field E can not only pump 
electrons but also accelerate them, as revealed in our 
gauge invariant kinetic Bloch equations [Eq. ©j^ 8 - by 
the pump and the drift terms, respectively. However, in 
previous works by Knorr et aZ. ) 33 ' 35 ~ — the drift term is 
not included and hence the case of low pump energy, for 
which the drift term is expected to be important, can not 
be investigated with their equations. In contrast, our ki- 
netic Bloch equations can deal with both the high and 
low pump energies completely. 

Under the laser field with pump-photon energy to l , the 
off-diagonal terms of the density matrices are generated 
around the resonant absorption point (where 2v F k ~ 
Wi). They oscillate with the frequency around ljl- Fol- 
lowing the approach adopted in semiconductor optics*^ 8 - 
we introduce the slowly varying interband-polarization 
component P M k = Piik.i.-ie 1 " 1 **. Then, the kinetic equa- 
tions can be rewritten into 

9t/ M k,„ = 77lm(0* k P Mk ) + |e|E • V k / M k,^ + 9t/ M k ,r/ (scat j 

dtP^ = -i£l(tk(fpk,i - //<k,-i)/2 - i^k-Pjuk 

+ |e|E- VkP M k + 9tP M k|scat, (10) 

with the detuning v£ = 2v F k — ujl + 
Tr K S "k + ^Pump)^] and the Rabi frequency 

n M k = -2e^ t (E™ , )lj _ 1 +^ nmP)lj _ 1 ) 

= p\e\v F E sm6 k e~^(l - e 2i » Lt )/w L - 2E^' )1) _ 1 e^ t . 

(11) 

When the pump-photon energy is high enough, these 
equations can be further simplified by neglecting 
the high-frequency oscillating terms (rotation-wave 
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approximation) and the Rabi frequency becomes 
= ^\e\v F E sm9 k e~ t2 ^ 2 ^/u L 

+ 2 E s k,k- q , 1>1 ^- q , k> -i,- 1 ^ r (q 5 o)p, k _ q . (12) 



In this work, the q-resolved dynamics of the RI and op- 
tical phonons are also studied, while the acoustic phonons 
are set to be in equilibrium with the environment at tem- 
perature To. The kinetic equations of the hot phonons 
are 



Then, one finds that, except for the drift and scattering 
terms, Eq. (fTOj) is of the same form as that in semicon- 
ductor optics4£ However, when the pump energy is low, 
the rotation-wave approximation is no longer valid and 
the Rabi frequency must be calculated with Eq. (ITT1) . 
We further define a new set of matrices p^ with 

PiJ.k,v,v = //ik,7; an d |3/ik,i,-i = Ptk,-i,i = ^k- Un- 
der such definition, p^ vary slowly with time and the 
scattering terms in Eq. (fTO)) . which include the electron- 
electron Coulomb, electron-impurity, electron-optical 
phonon and electron-remote-interfacial (RI) phonon 
scatterings, are given as 



0, 



t P fXk.,rjrj' 



E 

q»)i»)2 V3 



{[ n E ^k,k-q,7j7ji^k- 



q,i?i,»?2 



k— q,k,r/ 2 r]3> ^^ik, 7/3,77' 



E t/r (q,e r , ;kl+q - E^ki) 



+q,??4»7i'Vki+q,J7i,J7 2 



x V a (q, e^kx+q - £ v' 3 ^i)^k 1 ,k 1 

x S M ' d < 

ki +q,ki ,n' 2 lf 3 ^'ki ,r)' 3 ,r]' 4 
x e i(v-ril+m-ri3+'q' i -ri' 1 +'q! 2 -i-i! 3 )uj L t/2 

X 5(£^ 2 k-q — £?; 3 k + ^jjjki+q — e 773ki)] — 
+ {77 « > T]'}\ 



><]} 



(13) 



0, 



kf>flk,T]T]'\ e i — {[ 7 1 " E ^k.k-q.r/r/i^k-q,?)! 



'(2 



q»;i»?2 

X k-q,k,l, 2 l,2^k,t) 2 ,l'l^^l2 k £jj2k-qj| e 

x <5(e„ 2 k- q - £r, 2 k)] - [>^<]} + {v «— > V}*, (14) 



9* ft 



tk.,7]7)' 



{ 



E 

qm^ w'^± 



fifi' X 



Mi 



k,k — q,r/r/i k— q,k,r/2?7; : 



X (^q ±A ^k-q^^^k,, 3 ^ 



qPji'k-q.i]! ,?72 ftik,r; 3 ,?j' 



X e 



i(r]-rn+ij2-V 3 )uj L t/2 



+ {»7 <—>»/}*> 



(15) 



in which [>•< — ><] stands for the same term in the previ- 
ous [ ] but interchanging >< — X, and {r/ < — >■ 77'}* rep- 
resents the conjugate of the same term in the previous 
{ } except the exchange of 77 and 77'. In Eq. (fT5|). iV q A = 
n± q +l/2±l/2 with 7i q standing for the distribution func- 
tion of the phonon in branch A; V a (q, Ej/k+q — e r) k) = 
[U r (q,£,,' k +q-e^k)]*; P> k = l-p^k and /5< k = p^. The 
optical phonons include the transverse optical phonons 
(KO) near the K (K') point and the longitudinal (LO) 
as well as transverse optical (TO) phonons near the T 
point. The detailed forms of the scattering matrix ele- 



ments i/(£ r , 2 k - £r, 2 k- q ) and 
pendix B. 



are given in Ap- 



dm 



qlcp 



ft^qlpp, 



(16) 



in which d 1 



in 



qlcp 



come from the carrier-phonon scatter- 
ing and dtriq | pp describe the anharmonic decay of hot 
phonons in relaxation time approximation. They are 
given by 



qlcp 



k?)ir) 2 T;3774/i/i' 

X ( n > rt < N +X — 

A ^'k,?7i,t72^A'k+q,n3,'74 q 



k.k+q,?7 2 773 



fti'k,?)! ,772 Pfik+q 



mm' ' -q > 



x e 
qipp 



i(V2— m+Vi— Va)uiLt/2 



<5(£r )3 k+q - £r )2 k ~ W q \j\ , (17) 



pp' 



(18) 



where r pp is the phenomenological relaxation time from 
the phonon-phonon scattering and n° is the number of 
the A branch phonons at environmental temperature To. 

Before we show our numerical results, we first give a 
simple analysis on the Auger process. If only the diagonal 
terms of the density matrices (i.e., 77 = 77' , 771 = 772 , = 
rfi Vi — V2 and 773 = 774) are involved, the Coulomb 
scattering terms become 



'Aik- q ,77i 



q»?i 

x ^u.-q,k,mn^ k <v E U r (q, E^ki+q - £»^ki) 

X V (qj^ki+q — £t7 2 k 1 )^ki,ki+q,J7 2 Jj^/yki+q,Vi^//ki,J7 2 
X ^ki+q.ki^Jjj^^'Jik-q _ £ -')k + C^ki+q - £»; 2 ki)] 



- >< 



►< 



(19) 



From this equation, one finds that it describes the scat- 
tering process with one electron being scattered be- 
tween bands 77 and 771 while the other electron be- 
tween bands 772 and 77^ under the energy conservation 
condition £(E„ lk _ q - e„k + £^ kl +q - e^kj- More- 
over, the vertices of the involved two scattering pro- 
cesses, ^ r (q,^ lkl+q -£ J?2 ki)5'k 1 ,k 1 +q^ 2 T ; i' S, k,k-q, nni and 

V™(q,£ Vikl+q - ^) S L^, mv S k 1+ ^M,r,W 2 ' are con " 
jugated with each other. For the Auger-type scatter- 
ing process here, it means that three out of the four 
involved band indices (77, 771, 7/4 and rj' 2 ) are the same 
and the left one is different, corresponding to the case 
that one of the two scattered electrons is transferred 
from one band to the other while the other electron 
remains in the same band. Then, taking 771 = 77^ = 
V2 V f° r example and substituting the energy dis- 
persion, the requirement of energy conservation becomes 
v F (|k-q| + |k| + |ki+q|-|ki|) = 0. Since |k-q| + |k| > 
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|k-q-k| = q = |k x -h 
finds that the condition 

v F q 



ki| > |ki| - |ki +q|, one 



(20) 



must be satisfied. Although this condition is obtained 
in the special case 771 = r/[ = j]' 2 ^ 77, one can prove 
that it is valid as long as the band indices in the delta 
function is Auger type. On the other hand, it has been 
shown that, under this condition, the dielectric function 
under the dynamic screening e(q, e^ki+q — e ^ki) also 
diverges * 46 ' 47 Therefore, both V r (q, e^jki+q ~ e ?^ki) an d 
V a (c[, £, ; jk!+q — e r 1 ' 2 k 1 ) are reduced to zero and the cor- 
responding processes are forbidden. Nevertheless, if the 
scatterings involving the off-diagonal terms are taken into 
consideration, terms like 



TABLE I: Parameters used in the computation. 





a 


1.42 A 


Vf 


1 1 r\8 / a, 

lx 10 cm/s 




UJ-p 


lyo mev 


\ d t) 


40. OU ev /A 




UK 


161 meV b 


\ d k) 


92.05 ev /A 




D 


19 eV c 




2xl0 6 cm/s c 




pra 


7.6xl0~ 8 g/cm 2c 


d 


0.4 nm d ' e 


SiOa 




59 meV c 


91 


5.4 x 10" JC 




wri 2 


155 meV c 


92 


3.5 x 10" 2c 




r s 


0.8 f 


n Ic i 


1.5 s 


SiC 


ijJm 


116 meV h 


9 


1.4 x 10 _2h 




r s 


0.4 c4 




2.6 j 



' Refs. .2. b Refs. 70 and 71. c Refs. 72 and 73. d Refs. 62. 
'■ Ref. 63. f Refs. 51 and 61. g Ref. 74. 
Ref. 76. j Ref. 77. 



Ref. 75. 



X ^k-q,k,?727)3^Mk,i73,7)'^ (l) ^Jkx+q — e »i;ki) 

X V (q, £,^ki+q — ^akl^ki.ki+q.rj^'Vkl+q.l)! 



'/2 



■ / k 1 4-q,ki,)7^''/i'ki,j7^,r7 1 e 'fl'J'Ji 



only show the results in K valley (jx = 1) with the val- 
ley index /1 no longer appearing for simplicity. With the 
evolution of carrier distribution, the temporal evolution 
of the optical transmission at the probe-photon energy 



w pr is calculated fro m 13 ' 66 i 67 



T Wpt (f) = |1 + iViay(T Wpr (i)VMo/eo/(l + 



(22) 



Here, 8w 



indicates an 



also exist in c^p^k.jjjj' i cc - u r?i,->j;»)»)i 
Auger type scattering process for the scattering with 

the vertex V r (q, e^ki+q-e^kJSkV 



qij, 

-q,Wl ki,ki+q,77[?7j 

However, 5fj' 2 ^ V2r]3 guarantees that the scattering pro- 



cess restricted by the energy conversion [whose vertex 
is F a (q,£^ kl+q -£^ kl )5^ +q)kiiJ? ,-,^_ q!k)%7) J is non- 
Auger type, which means that the requirement Eq. (12TJ)) 
is no longer necessary. Therefore, these scattering 
terms are permitted and are expected to contribute 
to the dynamic of carriers. Moreover, for the corre- 
sponding scattering term [Eq. (|21[)]. there is a prefactor 
exp[?(?7 - 771 + fj[ - r}[)u L t/2] with 77 - 771 + fj[ - r)[ ^ 0. 
Since the other factors in the scattering term all vary 
slowly with time, the contribution of such scattering term 
oscillates very fast when wl is large and hence it has little 
effect on the dynamics of carriers. However, when lol is 
small enough, this kind of Auger process is expected to 
have marked influence on the evolution of electron distri- 
bution function. 



where n rc { is the refractive index of the substrate and 
iViay is the number of graphene layers. Here, the corre- 
sponding interband optical conductivity, which is deter- 
mined by carrier distributions, is given b y 13 i 54 i 68 i 69 
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d9^ sin 2 (0 k J[/k„a(i) - /k. 



in which |k w | = w pr /(2ui?) is the resonant absorption 
point. DT is calculated from AT Wpr (» /T^ = [T UpT (t) ~ 
T° ]/T° , with T® representing the transmission be- 
fore the pump. We present our numerical results with 
material parameters given in Table |H The full width 
at half maximum (FWHM) of the pump pulse is cho- 
sen to be 100 fs (corresponding to at — 42.5 fs) and 
T pp = 2.8 ps^ The equilibrium chemical potential is 
taken to be 18.6 meV (the corresponding electron den- 
sity is 1.43 x 10 11 cm -2 ), the environmental temperature 
is set to be 300 K and the impurities are absent, unless 
otherwise specified. 



III. RESULTS 



The kinetic Bloch equations, Eqs. ([6]) and (fl"6|) with 
all terms included are numerically solved following the 
method detailed in Ref. [60I except the drift term. In 
this paper, the drift term is treated up to the third or- 
der to suppress the noise (see Appendix C). Then, the 
temporal evolutions of the carrier and phonon distribu- 
tion functions are obtained. From Eq. (I10[) . one finds 
that / Mk ,r, = J- M k,r, and P^ k = -P-fik- Therefore, we 



A. High pump-photon energy 

We first investigate the dynamics of carriers and 
phonons with high pump-photon energy. The substrate 
is chosen to be Si02- In the calculation, ljl = 1500 meV 
(corresponding to the near- infrared light), unless other- 
wise specified. 
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FIG. 1: (Color online) (a)-(c) ln(f^\ — 1) as function of k 
with u)l = 150 meV along the x and y directions at different 
times. Here, the pump process is switched off. The black 
dotted curve in (c) is (i>p|k|— /4l)/(&bTi) with T\ = 333 K and 
Hi — 11.4 meV. To make the results clearer, the curve for t = 
100 fs plotted in (c) is ln(/^j — 1) + 1. (d) Temporal evolution 
of the total momentum k nct (t) with ujl = 150 meV (solid 
curve with the scale on the left hand side of the frame). The 
temporal evolution of the electric field —E(t) is also plotted 
for reference (dotted curve with the scale on the right hand 
side of the frame). 



1 . Influence of drift term 

We first investigate the influence of the drift term on 
the dynamics of electrons. To make it clearer, the pump 
process is switched off by excluding -ffp ump in Eq. © . In 
the calculation, the maximal electric field in Eq. (|3|) is 
taken as Eq = 150 kV/cm and the pump-photon energy 
ujl = 150 meV. The k dependences of \n(f^l — 1) along 
the k. x and k y directions at different times are plotted in 



Figs. [Ha) and (b), respectively. The dynamics of holes 
is similar and is not shown here. From the figures, one 
finds that before t — —100 fs the influence of the electric 
field is negligible. Moreover, the electron distribution at 
that moment is the Fermi distribution, with 



(7?v F |k| - fi ri )/(k B T n ) 



(24) 



as shown in the figure. Here, T v and fi^ are the elec- 
tron temperature and chemical potential in correspond- 
ing band rj, respectively. After t = —100 fs, the elec- 
trons arc drifted under the influence of the electric field 
[Eq. ([3])] along the x direction. From this oscillat- 
ing electric field, a net oscillating momentum k net (i) = 
Z)k k /k4(i)/Z)k/k,i(*) alon S tue x direction is gained 
by electrons. As a result, the location of the minimum 
of ln(/j~/J — 1) along the k^ directions oscillates around 
k = 0, as shown in Fig. [Ha). However, the minimum 
along the k y directions remains fixed at k = but its 
magnitude oscillates, as revealed in Fig. [IJb). To further 
illustrate the influence of the electric field on the net mo- 
mentum k no t(£), we plot the temporal evolution of k not 
in Fig. QJd). The electric field —E(t) is also plotted in 
the same figure for reference. One finds that its phase is 
always n/2 in advance of fc nct . Moreover, after the pulse 
(t > 100 fs, when \E\ < 0.1\E \), the net momentum 
gained from the electric field tends to zero. In addition, 
it is found that for curves at t = —100, and 100 fs in 
Fig. QJa), the magnitudes of slopes of the curves away 
from the minima decrease with the temporal evolution, 
indicating that electrons are heated by the drift pulse [see 
Eq. (|24p]. Then, with the negligible net momentum, the 
isotropic hot-electron Fermi distribution is expected to be 
established under the Coulomb scattering after the pulse. 
To show this, we plot the k dependences of ln(/j~/£ — 1) 
along the k x and k y directions at t = 100 and 130 fs 
in Fig. [ljc). One finds that at t = 100 fs, the curve is 
almost linear with k for each direction, which indicates 
the establishment of the Fermi distribution along that 
direction. To further show the buildup of the isotropic 
Fermi distribution, one fits the distribution with Eq. (|24|) 
along the direction 9^ and obtains the corresponding hot- 
electron temperature T v g k and chemical potential // r) e k . 
Then, the establishment of the isotropic Fermi distribu- 
tion is identified when T v g k along different directions are 
very close. Here, it is found that at t = 130 fs the rela- 
tive difference of the temperature increase T\g k — Tq along 
different directions is less than 10 %. 

Then, we investigate the influences of the pump- 
photon energy and impurity density on the thermaliza- 
tion of electrons by the drift term. We plot the pump- 
photon energy and impurity density dependences of the 
hot-electron temperature T\ in Fig. [2j With different 
ujl and Ni, the isotropic Fermi distributions are estab- 
lished at different times. Here, T\ is fitted along the k^ 
direction at t = 230 fs, when the isotropic Fermi dis- 
tribution has just been established (relative difference 
between temperature increase T\Q k — To along the k^ 
and kj, directions is less than 10 %] for the slowest case 



7 



Nj (10 12 cm" 2 ) 



10 



20 



30 



40 



50 



340 



330 



320 



310 



300 



(Dl = 300 meV — >- - 

500 meV 

Xjc" 


t = 300 k ; 

t = 230 fs - 






ir , 



480 
450 
420 
390 
360 
330 
300 



150 



250 350 
C0 L (meV) 



450 



FIG. 2: (Color online) Pump-photon energy lul and impurity 
density Ni dependences of the hot-electron temperature T\. 
Here, the hot-electron temperatures are obtained at t = 230 fs 
and the environment temperature is To = 300 K. The impu- 
rity density dependence of Ti are calculated with oil — 300 
and 500 meV and its scales are on the top and right-hand side 
of the frame. 



in our calculation. The impurity density dependence 
of Ti is obtained with wl = 300 and 500 meV. It is 
shown that T\ decreases monotonically with the increase 
of the pump-photon energy. Moreover, a peak appears 
in the impurity density dependence as shown in Fig. [2] 
These phenomena can be well understood from the Boltz- 
mann equation under the relaxation time approximation: 

d t fk,r, = |e|E ■ Vk/k,^ - (/k,„ - /k,^)/ T p ; with fk,n bein S 
the equilibrium distribution and t p standing for the relax- 
ation time. From this equation, one obtains the conduc- 
tivity under the AC electric field E = Eq exp (— iuiLt)x 
as 



2, ,2 



e v 



dkk- 



1 + (uiLTp) 



2 £( 

V 



dfk, v 



de. 



77k 



, (25) 



which is the same as the intraband conductivity given in 
the literature! 39 ' 54 ' 67 From this equation, it is found that 
the conductivity and hence the energy absorbed from the 
electric field decreases monotonically with the increase 
of uil and a peak exists in its N (or equivalently 1/t p ) 
dependence. These phenomena can also be understood 
more physically. One finds from Eq. (|2"5j) that the con- 
ductivity tends to zero in the limit of no scattering or 
infinitely large scattering. This can be comprehended 
as the effect of the electric field on the electron is to- 
tally cancelled out in one oscillation period in the case 
without scattering and the electrons can not be drifted 
by the electric field in the limit of infinitely large scat- 
tering strength. Then, in the case with finite scattering 
strength, a peak must appear as shown in Fig. [2J As for 
the decrease of T\ with the increasing w^, it is under- 



stood as the electron becomes more and more difficult to 
follow the laser field with the increase of the oscillation 
frequency. 

To further compare the influence from the drift and 
the pump terms, we also perform the calculation with 
the pump process included but the drift term excluded. 
In the calculation, u>l = 500 meV, N — and the other 
parameters remain the same. At the time in Fig. [2] i.e., 
t = 230 fs, the isotropic Fermi distribution is checked 
to have already been established and the obtained hot- 
electron temperature is Ti = 854 K, much larger than 7 K 
established by the drift term only. In the sample usually 
reported in the experiments ; 79 ' 80 the mobility is larger 
than 1000 cm 2 /Vs {N < 2.5 x 10 12 cm~ 2 ). Within this 
impurity density range, the temperature increase due to 
the drift term is less than 37 K. This is much smaller than 
that from the pump process (the influence of impurity on 
the pump process is small with such pump intensity and 
impurity density range). Moreover, since the thermaliza- 
tion due to the drift effect decreases with increasing lol, 
the drift term is negligible in the case with pump-photon 
energy higher than 500 meV. This is consistent with the 
investigation in semiconductors where the energy of the 
pump photon is very high due to the large band gap4^ 
However, in the case with very low pump-photon energy 
as we will investigate in the next subsection, the drift 
term is expected to be very important, especially for the 
case where the electric field does not change direction 
during the pulse and a large net momentum is transferred 
to the electrons via the drift term. 



2. Dynamics of electrons 

Before the investigation on electron dynamics, we first 
calculate the DT and compare it with the experimental 
data in single-layer graphene by Hale et aL 25 [Fig. [3]Ja)]. 
It is noted that this experiment has been fitted in our 
previous worki 3 ^. However, in that paper, the Fermi dis- 
tribution was assumed to be established at all the time 
and the hot phonons are averaged with a fitting param- 
eter -E ma x7 which describes the upper energy limit of the 
hot carriers being able to emit phonons. In this paper, 
with the inclusion of the pump term and the q-resolved 
hot phonons, the buildup of the Fermi distribution can be 
studied elaborately and the fitting parameter E max is no 
longer needed. Then, the influence of the approximations 
in our previous paper can be checked. In the calcula- 
tion here, the probe-photon energy u; pr = 1100 meV, the 
FWHM is taken to be 180 fs and the pump-photon en- 
ergy lot, — 1500 meV, as indicated in the experiment. It 
is noted that with N = and such high pump-photon en- 
ergy, the drift effect of the laser field is negligible as shown 
in the previous subsection. In addition to these parame- 
ters, the fitting parameter t pp = 2.8 ps and the maximal 
electric field Eq — 300 kV/cm which leads to the maxi- 
mal electron density being 2.7 x 10 12 cm~ 2 . This electron 
density is comparable with the experimental estimation 
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FIG. 3: (Color online) (a) DT from the numerical results com- 
pared with the experimental data in single-layer graphene. 25 
The results are normalized as Hale et al.r^ The arrow in the 
figure indicates t — 133 fs, which is the beginning of the re- 
laxation. (b)-(d) ln(/,7j — 1) as function of k along the x and 
y directions at t — —90, 0, 100 and 131 fs. The black dotted 
curves in (b) and (c) are the case before the pump pulse. The 
curve for t = 100 fs plotted in (d) is ln(fe{ - 1) + 1. The 
cases without the Auger process are also plotted in (b)-(d). 



and that in our previous paper (4.6 x 10 12 cm -2 ). For 
the relaxation of DT which begins from t = 133 fs, as 
shown in Fig. |3{a) , one finds that our numerical results 
agree very well with the experimental data with a fast re- 
laxation followed by a slow one with the relaxation times 
being 0.28 and 1.33 ps, respectively. The discrepancy be- 
tween our numerical results with the experimental data 
at the initial increase of the DT is understood to come 
from the approximation in our model where the inter- 
ference between the pump and probe pulses is neglected 
and the Markovian approximation^ is applied. 

Then, we investigate the dynamics of electrons, espe- 
cially the buildup of the Fermi distribution under the 
linearly polarized light. We plot in.(f^\ — 1) as func- 



tion of k along different directions at different times in 
Figs. OHb)-(d). With the strong scattering and the neg- 
ligible drift effect, the distribution is center symmetric 
and hence the distribution functions along the — x and 
— y directions are not shown. Besides, the dynamics 
of holes is almost the same^I and is not plotted here. 
From Fig. G^b), one finds that the distribution of the 
pumped electrons are anisotropic and mainly around 
k = 1.14 nm -1 ss uil/(2?ivf) along the y-axis as shown 
by the valley there. This comes from the linearly polar- 
ized pump light as shown by sin (9k in Eq. (jlip . During 
the pump process, the scattering also plays an important 
role by smearing out the valleys and making the distri- 
bution tend to be isotropic as shown in Figs. [3Jb) and 
(c). Moreover, it is found that the slopes of the curves 
around the Dirac point decrease with the temporal evo- 
lution, indicating their thermalizations under the scat- 
tering. At t = 131 fs shown in Fig. [3Jd), one finds that 
the curves become almost linear with k for each direc- 
tion. By fitting the curves with Eq. (|24p. the obtained 
hot-electron temperature Tig k along the x and y direc- 
tions are 2160 and 2221 K, respectively with the relative 
difference for these two directions being less than 3 %. 
By noticing that the pulse [Eq. d3j] ends after t = 165 fs 
(when \E\ < 0.1|£o|), one concludes that in the case with 
such long pump pulse width, the hot-electron Fermi dis- 
tribution is established during the pulse and once it is 
established, it is an isotropic one. Considering that the 
relaxation of the DT begins at 133 fs which is after the 
establishment of the Fermi distribution, the approxima- 
tion of the buildup of the Fermi distribution in our pre- 
vious work is acceptable. However, this approximation 
leads to the inaccuracy of the hot-electron temperature 
(3200 K at the beginning of the relaxation in our previous 
work—) . 

To further show the influence of the Auger process in- 
volving the interband coherence, we also plot the results 
without the Auger process in Figs. G3b)-(d). One finds 
that the exclusion of the Auger process has little influ- 
ence on the evolution of the distribution function. This is 
in agreement with our argument in Sec. II that the Auger 
process which survives the dynamic screening in the pres- 
ence of the optical polarization does not contribute due 
to the high pump-photon energy. Similarly, also due to 
the high pump energy, the rotation-wave approximation 
is checked to be valid (not shown in the figure). 



3. Dynamics of phonons 

With the inclusion of the q-resolved hot phonons, we 
are able to give a more detailed investigation on the dy- 
namics of phonons than our previous work^i We calcu- 
late the temperature T p h(q) for phonons with momen- 
tum q at branch A by T p h(q) = oj\/[kB ln(l + 1 /n^)] and 
then, plot them in momentum space at t = 131 fs for 
the KO and Mi phonons in Figs. Ufa) and (b). The q- 
resolved temperatures of the other phonons are similar to 
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FIG. 4: (Color online) (a) and (b) Typical q-resolved hot-phonon temperatures of KO and RIi phonons at t = 131 fs, 
respectively. The black dotted circles indicate q — oj q \/vF- (c)-(f) The angle averaged hot-phonon temperatures for different 
phonon branches as function of VFq at t = —90, 131, 300 and 600 fs, respectively. The arrows indicate the location of the 
corresponding phonon energy VFq = u) q \. It is noted that the phonon energies of TO and LO phonons are degenerate and we 
only plot the arrow for TO phonons here. The region close to VFq — is enlarged in the inset for KO and RI2 phonons at 
t = 600 fs. 



these two phonon branches. It is shown that due to the 
anisotropic distribution of electrons, the thermalization 
of phonons under the electron-phonon scattering is also 
anisotropic. Moreover, it is found that the thermaliza- 
tion of phonons with small q (q < uj q \/vp, which are the 
phonons within the black dotted circle) is negligible. 

This can be understood from the energy conservation 
during the electron-phonon scattering process [Eq. (JTSJ)], 
J(£, 2 k- q -£i, 3 k±w 5 A)- By using the relation vp\\t — q| — 
up|k| < VF|k — q — k| = vpq, one finds that the phonons 
with q < uj q \/vF and q > uj q \/vp are involved in the 
interband (772 7^ 773) and intraband (772 = 773) electron- 
phonon scattering, separately. Therefore, the negligible 
thermalization of the phonons with q < uj q \/vp means 
that before the establishment of the hot-electron Fermi 
distribution (t < 131 fs), the recombination of electrons 
and holes due to the interband electron-phonon scatter- 
ing is marginal. In order to show this clearer, we plot 
the q dependence of the angle averaged phonon tem- 
perature TP h = J Q 2Tr d#qT p h(q)/(27r) at different times 
in Figs. EJc)-(f), as done by Knorr et al.M^ The loca- 
tions vpq — u q \ for each branch of phonons are indicated 
by arrows in Figs. Htc)-(f). It is found that in the first 
several hundred femtoseconds, only the thermalization of 



the phonons with vpq > u q \ is efficient. After that, the 
temperatures of the phonons with vpq < u q \ increase 
markedly. Especially, for the KO and RI2 branches, 
the hottest phonons appear among the phonons with 
vpq < u q \ as shown in the inset of Fig. IDJf). This is 
consistent with the results by Knorr et a L 32 ^ 35 and can 
be understood as the interband electron-phonon scatter- 
ing happens only among the low energy electrons with its 
energy |e^k| < u q \. However, the energy of the photoex- 
cited electrons are very high {lul/2 > 3.9uj q \) and before 
the low-energy electrons are heated under the scattering, 
the energy gained by the electrons first relaxes through 
the intraband electron-phonon scattering process. Then, 
after the thermalization of the low energy electrons, the 
heating of the phonons with VFq < u q \ starts to be im- 
portant. 

One also observes from Figs. Ufa) and (b) that the 
hottest phonons are along the q y direction. This comes 
from the angular dependence of the electron-phonon scat- 
tering matrices. Since the hottest phonon appears with 
q > u q \/vF, we only consider the intraband scatter- 
ing process. Then, with only the diagonal terms of 
the density matrices taken into consideration, the scat- 
tering matrices in Eq. (fllpl) l-^k k- q w l 2 are P ro P or - 
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tional to 1 — cos(# k — #k- q ) for the KO phonons and 
l + cos((9 k — #k-q) for the RIi phononSi 63 i 70 ' 71 One imme- 
diately understands that the parallel or antiparallel scat- 
tering are the strongest. Then, since the photoexcited 
electrons are mainly along the k a direction, the hottest 
phonons, which come from the elcctron-phonon scatter- 
ing, also have to be along the y direction as shown in the 
figures. In addition, Figs. Hfc) -(f) show that the thermal- 
ized phonons in KO and LO branches occupy larger mo- 
mentum range than the other phonons. This also comes 
from their different angular dependences of the scattering 
matrices as revealed by Malic et al.£5- 



B. Low pump-photon energy 

We then investigate the dynamics of electrons and 
phonons with the pump-photon energy down to the 
THz regime. The maximal pump field in Eq. §3§ is taken 
as E = 5 kV/cm. The substrate is again chosen to be 
Si0 2 . 



1. Dynamics of electrons 

In this subsection, we investigate the dynamics of elec- 
trons with lul — 15 meV (3.6 THz). With such low 
pump-photon energy, the corresponding oscillation pe- 
riod of the electric field is 276 fs, which is even larger 
than twice of the FWHM. Therefore, during the pump 
process, the laser field E almost keeps along one direction 
and a large net momentum is transferred to electrons. In 
this case, the drift term is very important since it de- 
scribes the momentum transfer. Furthermore, due to the 
transferred net momentum, a drifted Fermi distribution 
after the pulse is expected. In previous works with static 
electric field) 45 ' 81 it has been shown that in the case with 
strong Coulomb scattering, the drifted Fermi distribution 
is expressed as 

h n = {exp[(?7« F fc - u„ • k - ^)/{k B T ri )] + l}" 1 , (26) 

with u n being the drift velocity. However, in graphene, 
it has been shown by Zhou and Wu that the Coulomb 
scattering is not strong enough compared with the drift 
effect of the static electric field and the drifted Fermi 
distribution is instead described by^ 5 - 



/ k?) = {exp[(?7u F |k - u n | - (I,,) I \k B T v )} + 1}~ 



(27) 



Although both distributions are the same in semiconduc- 
tors with parabolic spectrum, they have different prop- 
erties in graphene due to the linear dispersion. For 
Eq. ((271), one finds that u, ; = S k k/ k>J7 (t)/^k/k,»?(*) is 
nothing but the center-of-mass drift velocity. However, 
in Eq. (|21)1) . is only an effective parameter. Moreover, 
in the k dependence of ln(/j7} — 1) along the direction 
of u^, Eq. (|26[) indicates a minimum at k — and the 
magnitudes of the slopes for # k and # k + tt directions are 
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FIG. 5: (Color online) (a) V^if^l, ~ 1) along # k = and 7r 
directions for conduction band electrons (jj = 1) at t = —50, 
25, 100, 300 and 500 fs, as well as that for valence band 
electrons (r) = —1) at t = 300 fs. Here, ujl = 15 meV and the 
impurity density Ni = 0. The black dotted curve is calculated 
from Eq. (26j) with Xi = 323.6 K, ui = 0.12^ F and = 
16.3 meV. (b) Temporal evolution of T v e k and u v e k (with 
the scale on the right-hand side of the frame) for electrons in 
conduction (solid curves) and valence (dashed curves) bands 
obtained by fitting Eq. ()26l) along directions 6k — and n 
(red squares), as well as # k = 7r/4 and 57r/4 (blue dots). 



different. However, Eq. (1271) shows that the minimum 
deviates from k = and the magnitudes of the slopes are 
identical. To reveal the distribution established here and 
the dynamics of electrons, we plot h^/^ — 1) as func- 
tions of k along the direction of the electric field E (i.e., 
9 k = and ir) at different times in Fig. EJa). One finds 
that during the pump process, the minimum of the curve 
is first drifted away from the Dirac point as shown by the 
curves at t = —50 and 25 fs. However, before the pulse 
ends (t > 100 fs, when \E\ < 0.1|-Eo|)j the minimum has 
already started to tend to the Dirac point and the magni- 
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tudes of the slopes for the curves away from Dirac point 
become different. Then, at t = 300 fs, the minimum is 
around the Dirac point and the curves along #k = and 
7r directions are almost linear with k, but with different 
magnitudes of slopes. With such different magnitudes of 
slopes, the corresponding drift of the center-of-mass fc n ct 
is 0.02 Dm" 1 . This indicates the establishment of the 
distribution Eq. (1261) . as shown by the fitting plotted as 
the black dotted lines in the figure. This is understood 
that the drifted Fermi distribution is established after 
the pulse in our present investigation. Then, with the 
Coulomb scattering being the dominant scattering, the 
distribution as Eq. (l26l) is established. 

It is also shown in Fig.JSfa) that at t = 300 fs, the sim- 
ilar drifted hot-electron Fermi distribution in the valence 
band (r) = —1) has also been established since the mini- 
mum of the corresponding curve is around k — and the 
curves along 9^ = and ir directions are almost linear 
with k with different magnitudes of slopes. Moreover, 
the distributions along the other directions #k are also 
similar. To further investigate the buildup of the uni- 
tary drifted Fermi distribution along all directions, we fit 
the distributions of electrons along different directions in 
both the conduction and the valence bands with Eq. ([2"B| . 
The obtained T v g k and u v g k are plotted in Fig.JSfb). One 
finds that after t = 300 (380) fs, the relative difference 
between the drift velocities along different directions for 
electrons in conduction (valence) band is less than 5 % 
and the relative difference of the temperature increases 
T Wk ~ T (T_ie k - T ) is less than 10 %. This indicates 
the buildup of the unitary drifted Fermi distribution as 
Eq. (|26[) for all directions in conduction (valence) band. 
It is noted that the case studied here is n-doped, which 
leads to different buildup times, hot-electron tempera- 
tures and drift velocities of the unitary drifted Fermi dis- 
tributions for electrons in conduction and valence bands. 
However, under the interband Coulomb scattering, the 
temperatures and drift velocities in conduction and va- 
lence bands tend to be identical. From our calculation, 
the relative difference of the drift velocities is less than 
10 % after t = 650 fs and that of the temperature in- 
creases is less than 10 % after t = 350 fs, as shown in the 
figure. 

We also investigate the influence of impurities and 
pump-photon energy on the drift velocity of electrons. 
uie k obtained from the distributions along 0^ = and 
7r directions under different pump-photon energies and 
impurity densities are plotted in Fig. |U(a). It is noted 
that after t = 300 fs, the relative differences of the drift 
velocities along different directions are less than 5 % 
for all cases investigated here. One finds that in the 
case without impurity, the drift velocity Uig k and hence 
the total momentum relax slowly (with the relaxation 
time being about 1.8 ps). However, when the impurities 
with Ni = 5 x 10 11 cm~ 2 (the corresponding mobility 
is about 5000 cm 2 /Vs) are introduced into the system, 
the relaxation of the total momentum becomes very fast 
(with the relaxation time being 46 fs) and at t = 300 fs, 
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FIG. 6: (Color online) (a) Temporal evolution of «ig k for the 
cases with oil = 15 and 30 meV without impurity, as well 
as cj£ = 15 meV with impurity density N t — 5 x 10 11 cm~ 2 . 
(b) Temporal evolution of the photoexcited electron densities 
for the cases with all the processes included, without Auger 
process, without pump process and under the rotation-wave 
approximation (RWA), separately. Here oil = 15 meV and 
Ni = 0. 



the corresponding u%s k — O.OOSvp is negligible. This is 
because that the relaxation of the drift velocity comes 
from the electron-impurity and electron-phonon scatter- 
ings. However, the electron-phonon scattering strength 
is weak and only when the impurity density is large (e.g., 
5 x 10 11 cm -2 ), can the relaxation of the net momentum 
become fast. Besides, it is also shown in the figure that 
when Wi increases from 15 to 30 meV, the value of Uig k at 
t = 300 fs decreases markedly. This is because that with 
ujl = 30 meV, the laser field E changes direction during 
the pulse and the net momentum gained by electrons is 
partially cancelled out. 

We further investigate the influence of the pump and 
Auger processes on the evolution of the photoexcited 
electron density by switching them off separately. The 
results are plotted in Fig. [Hb). In the calculation, 
ujl = 15 meV and Ni=0. It is noted that without the 
pump process, electrons can still be excited from the va- 
lence band to conduction band. This is because that 
the electrons are heated by the drift effect of the laser 
field, which leads to the interband electron-phonon scat- 
tering. However, one finds from the figure that the dom- 
inant mechanism for the excitation of the electrons is the 
pump process in the case investigated here. Moreover, 
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it is also found that the Auger process indeed affects the 
dynamics of electrons and it leads to the change of the ex- 
cited electron density at the end of the pulse (t — 100 fs) 
as large as 22 %, which is consistent with our analysis 
in Sec. II. Furthermore, with such low pump-photon en- 
ergy, the widely accepted rotation- wave approximation in 
semiconductors optics is no longer valid. The temporal 
evolution of the photoexcited electron density under the 
rotation wave approximation is also plotted in Fig. [6jb) . 
One finds that the fast oscillation of the electron den- 
sity during the pump process disappears and the excited 
electron density at the end of the pulse changes about 
17 %. 



2. Dynamics of phonons 

We further investigate the dynamics of phonons. We 
plot the q-resolved temperatures of the hot RIi and KO 
phonons with wl = 15 mcV under the impurity densities 
Ni = and 5 x 10 11 cm -2 in Fig. [7] The other phonons 
are similar and not shown. The temperatures of these 
two branches of phonons are obtained at t = 300 fs when 
the unitary distribution of electrons for all directions has 
been established. From the figures, one finds that in the 
case without impurity, the q-resolved temperatures of the 
hot phonons are highly anisotropic. Moreover, by consid- 
ering that u,, in both conduction and valence bands are 
along the 9^ = direction, it is shown that the hottest 
phonons appear in the direction parallel to and the 
coldest phonons antiparallel to u v . This phenomenon is 
also shown in semiconductors under static electric field. 82 
It can be understood from the generation rate of phonons 
due to the electron-phonon scattering [Eq. ([T7|l] after the 
establishment of the drifted Fermi distribution of elec- 
trons. With only the diagonal terms of the density ma- 
trices taken into consideration, Eq. (|17|) is given by 



cWI =27rV 

vllcp / J 



iM,r i2j 



k+q,k,r; 2 »?i I ^( e »?2k+q £))ik w <?a) 
X (//T'k, t;i / A fk+q,r ;2 ^ V q FA - /^'k,» h /^k+q,l) 2 ^-q ) • ( 2 ^) 



Substituting the 
Eq. 126) and the 



drifted Fermi 
Bose distribution 



distribution 
of phonons 



riq = l/{exp[w g x/(^sJph)] — 1} into it, one finds 
that whether the phonons are emitted or absorbed is 
determined by the sign of the term 

e (j7iiipfe— U'k- fi+uiqi) I \k B T) _ e [j72t)p|k+q|— u-(k+q) — /i]/(fe B T) 

(29) 

Here, to simplify the analysis, the temperatures of 
phonons and electrons are taken to be the same T 
and the chemical potentials and drift velocities for elec- 
trons in different bands are also set to be unified fi 
and u, respectively. Then, due to the energy conser- 
vation <5(e^ 2 k+ q — e^ik — u q \), Eq. (|29|) is simplified 

to e (')i'i>Ffc-u-k-^+w,j A )/(/c£iT)j-L _ e -u-q/(ft fl T)l Q ne j m _ 

mediately understands that the hottest phonons appear 
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FIG. 7: (Color online) (a) and (b) [(c) and (d)] q-resolved 
hot-phonon temperatures with Ni = (5 x 10 11 ) cm -2 at 
t = 300 fs for RIi and KO phonons, respectively. 
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among phonons with u parallel to q and the coldest 
phonons appear among the phonons with u antiparallel 
to q as shown in Figs. [Tfa) and (b). This phenomenon 
can also be understood physically from the relaxation of 
the drift velocity of electrons. During the relaxation, the 
net momentum gained from the laser field is transferred 
from the electrons to the phonons, which implies the ten- 
dency of the absorption of the phonons with q • < 
and the emission of the phonons with q ■ > 0. On 
the other hand, in the system where the net momentum 
relaxes very fast through the electron-impurity scatter- 
ing, the distributions of phonons become more isotropic 
as shown in Figs.^c) and (d). 



C. Negative DT in system with high Fermi energy 

In this subsection, we further investigate the negative 
DT due to the weakening of the Pauli blocking. Sun 
et at measured the DT with high probe-photon energy 
(> 500 meV) and showed that a peak exists in the tem- 
poral evolution of DTJ^ Moreover, it is also found that 
there is a range of the probe-photon energy where the 
negative DT appears during the evolution. To under- 
stand this, they suggested that their epitaxial sample 
consists of a stack of graphene layer s 13 ' 83 ' 84 with some 
layers heavily doped and some layers undoped. More- 
over, the Fermi energy of the heavily doped layers is as- 
sumed to be around half of the upper boundary of the 
probe energy range for the appearance of the negative 
DT. Then, with the energies of the resonant absorption 
states of the probe pulse smaller than the Fermi energy 
in the heavily doped layers, the corresponding distribu- 
tion difference /k^.i — /k^.-i in Eq. (|2"B1 increases in the 
undoped layers but decreases in the heavily doped lay- 
ers after the pulse. This leads to the decrease of the 
conductivity in undoped layers but increase in heavily 
doped ones. These two opposite tendencies compete in 
the evolution and are responsible for the peak and the 
negative DT. Moreover, with the decrease of the probe- 
photon energy, the variation of the conductivity from the 
undoped layers is strengthened whereas that from the 
heavily doped layers is weakened. Furthermore, if the 
probe-photon energy is higher than twice of the Fermi 
energy in the heavily doped layers, the conductivities in 
both the heavily and undoped layers decrease. These 
properties lead to the fact that the negative DT appears 
only in a probe-energy range smaller than twice of the 
Fermi energy. 

To check above conjecture, we fit their experimen- 
tal data with our microscopic approach. The substrate 
is chosen to be SiC, the environmental temperature is 
50 K, the FWHM of the pulse is 100 fs and the pump- 
photon energy = 1550 meV as indicated by Sun et 
al.r^ The maximal electric field in Eq. <j3j> is taken to 
be Eq = 175 kV/cm. With such pump intensity, the 
photoexcited electron density is about 4 x 10 11 cm~ 2 
per layer. The probe-photon energy w pr = 550 meV 



is the same as that in the experiment. It is noted 
that with such high probe energy, the intraband absorp- 
tion is suppressed. 68 To model the multilayer structure 
in their sample, A^ayO^ (i) in Eq. (|22l) is replaced by 
Nh(Th(t) + N u a u {t) with au{t) and a u (t) being the con- 
ductivity in heavily and undoped graphenes, respectively 
and Nh and N u being the corresponding layer numbers. 
In our computation, ah (t) and a u (t) are calculated sepa- 
rately with the corresponding equilibrium Fermi energy 
being fi — 350 and meV, respectively, as estimated in 
the experiment. Nh = 2 and N u — 20 are taken as fitting 
parameters which are very close to those estimated in the 
experiment. The results are plotted in Fig. Ufa). One 
finds that our numerical results repeat the peak and the 
negative DT shown by the experimental data. Especially, 
when t > 500 fs, the fitting is quite well. The discrepancy 
between our numerical results with the experimental data 
for the magnitude of the peak is again from the approx- 
imation mentioned in Sec. IIIA2 as well as the simplifi- 
cation of the complex layer configuration to two kinds of 
doped layers. The electron distributions before and after 
the pulse are plotted in the inset of Fig. Ufa), which is 
consistent with the arguments by Sun et al.J^ Moreover, 
when the probe-photon energy is small enough or larger 
than twice of the Fermi energy, the DT keeps positive 
as shown by the curves with w pr = 520 and 710 meV in 
Fig. Elh). Based on these results, our microscopic inves- 
tigation supports their argument that the negative DT 
comes from the weakening of the Pauli blocking in the 
heavily doped layers. 



IV. SUMMARY 

In summary, we have investigated the nonequilibrium 
dynamics of carriers and phonons in graphene under the 
linearly polarized near-infrared and THz laser pulses via 
the microscopic kinetic Bloch equation approach i 48 ' 49 
The carrier-impurity, carrier-phonon and carrier-carrier 
Coulomb scatterings are explicitly included and the dy- 
namic screening for the Coulomb scattering is utilized. 
Moreover, based on the vector potential gauge and the 
gauge invariant approach^ both the drift and pump 
terms are naturally included in our approach. 

We first investigate the thermalization of electrons due 
to the drift term. It is shown that the thermalization is 
weakened with the increase of the pump-photon energy 
and it is even negligible when the pump-photon energy 
is higher than 500 meV. Moreover, a peak appears in 
the impurity density dependence of the thermalization. 
Then, the dynamics of carriers and phonons with the 
pump-photon energy up to the near-infrared regime is 
studied. We show that the temporal evolution of the DT 
reported by Hale et al~. can be well fitted. Moreover, 
under the linearly polarized laser pulse, the electrons are 
photoexcited anisotropically. However, their distribution 
tends to be isotropic under the scattering. As a result, 
the isotropic hot-electron Fermi distribution is found to 
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FIG. 8: (Color online) (a) DT from the numerical results com- 
pared with the experimental data extracted from Fig. 3(a) in 
Ref. The inset shows the electron distributions before 
(t = —100 fs) and after [t — 200 fs) the pulse in the heavily 
doped and undoped layers, (b) Temporal evolution of DT 
under different probe-photon energies. 



be established within t = 131 fs, which is before the end 
of the pulse (the FWHM of the pulse is 180 fs). This 
is consistent with the previous works i 26 ' 31 ' 35 In contrast 
to the high pump-photon energy case, in the case with 
the pump-photon energy down to the terahertz regime, 
we show that a large net momentum is transferred to 
electrons through the drift term and a drifted Fermi dis- 
tribution different from the one under static field is es- 
tablished in several hundred femtoseconds. 45 We further 
show that the Auger process investigated in the previ- 
ous work a 33 ' 34 i 36 which only involves the diagonal terms 
of density matrices is forbidden by the dynamic screen- 
ing. However, the Auger process involving the interband 
coherence contributes to the dynamics of carriers. Nev- 
ertheless, its influence is important only when the pump- 
photon energy is low. In contrast, for the rotation-wave 
approximation which is widely accepted in semiconduc- 
tor optics, we show that it fails when the pump energy is 
low. We also apply our microscopic theory to investigate 
the negative DT in the case where the equilibrium Fermi 
energy is high and the probe-photon energy is around 
twice of the Fermi energy. We show that the negative 
DT appears due to the weakening of the Pauli blocking 
under the pump pulse, which supports the suggestion by 



Sun et al. in their experimental workji 4 - 

In addition, the momentum-resolved hot-phonon tem- 
peratures are also studied. When the pump-photon en- 
ergy is high, the hot phonon temperatures are anisotropic 
due to the linearly polarized laser field and the hottest 
phonons appear with their momentum in the direction 
perpendicular to the laser field E. Moreover, we also 
show that the electron-hole recombination due to the 
interband electron-phonon scattering is marginal in the 
first several hundred femtoseconds. On the other hand, 
when the pump-photon energy is low and the relaxation 
of the net momentum is slow, the hot-phonon temper- 
atures are also highly anisotropic with the hottest and 
the coldest phonons appearing in the directions parallel 
and antiparallcl to the drift velocity, respectively, due to 
the drifted Fermi distribution of electrons. However, if 
the relaxation of the net momentum is fast as in the case 
with high impurity density, the hot-phonon temperatures 
become more isotropic. 
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Appendix A: Gauge Invariant Kinetic Bloch 
Equations 

We use the nonequilibrium Green function approach to 
derive the gauge invariant kinetic Bloch equations under 
the optical fieldi^^ The effective Hamiltonian without 
external field in the real space is given by^4 

C(-«V) = V F {~i l XT X V X - ITyVy), (Al) 

and the electron-optical-field interaction Hamiltonian is 
described byS 4 - 

H%£L,M = \e\v F [^r x A x (t) + r y A y (t)]/c. (A2) 

Then, in the base set of the eigenstates of Hq t (—iV), the 
above Hamiltonians read 

h o,t(p) = VFP&Z, (A3) 
#photo„,r(P^) = \e\v F A(t) ■ [pa z + M z x pa v ]/(cp),(AA) 



with p x = -id x , p y = -id y and p = ^Jp x + p 2 y . It is 

noted that this electron-optical-field interaction Hamil- 
tonian is consistent with the previous works based on the 
tight-binding approach j 55 i 56 
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The contour-ordered Green function 0^(1,2) are ma- 
trices with the matrix elements defined as G^ ,(1,2) = 

-i{T c [^{l)^{2)]). Here ^ and iptf are the Fermion 
field operators in the Heisenberg picture for electrons in 
band r] and valley fi; (1) = (ri,ti) and (2) = (r 2l t 2 ) 
stand for the space-time points on the contour and Gq 
represent the free-particle Green function. By multi- 
plying id ty - Ht} tT ($i) and -id t2 - H^ t (-^p 2 ) to 
the Dyson equation of G M (1,2) and subtracting the ob- 
tained two equations from each other, one obtains the 
equation with rj, ti, r 2 and t 2 . Then, transforming 
the equation to the center-of-mass and relative variables 
R = (r x + r 2 )/2, T = (h+t 2 )/2, r = n - r 2 and 
t = t\ — t 2: and omitting R due to the uniformity of the 
system studied here, one obtains^ 

id T G»(r, T, r) = [ff£ r (p r ) + ^ hoton , r ( P ,, T + r/2)] 

xG"(T ) T ) r)-G"(r I T ) r)[ff ,r(tr) 

+ H^ oton ^ r ,T-r/2)]. (A5) 

Here p r represents the momentum operator with respect 
to the relative variable. The gauge invariant functions 
are constructed as^ 



sound velocity and the corresponding scattering matrices 
are 



M k%\kr, = D Vq S k>ryMri S V.'li/\/ 2 P™ V Ph, (Bl) 

in which D is the deformation potential and p m denotes 
the graphene mass density 72 ' 73 For the RI phonons, 



= sf9VFe-« d Su v ,^5^/Jaq^, (B2) 



where g represents the dimensionless coupling parameter 
depending on the material of the substrate , 63 ' 75 a is the 
C-C bond distance and e^ 1 is the substrate-dependent 
dielectric function. For the two branches of phonons RT 
and RI 2 in the sample on Si0 2 , the corresponding dielec- 
tric functions are^ 5 - 

l/e™ 1 =l/[e i + e(q,0)]-l/[e o + e(q,0)], (B3) 
l/e™ 2 = l/[e oo + e(q,0)]-l/[e i + e(q,0)], (B4) 

with ei = 3.36, e = 3.90 and eoc = 2A0M For the only 
branch of RI phonons in graphene on SiC substrate, it 

:„85 



G"(k,T,T)= / ^rexp(-ik-r)G Al (r,T,r) 



(A6) 



with k = k /^ 2 dA|e|A(T + Ar)/c. Applying this 
transformation to Eq. (|A5I) and setting r = 0~, one has 



0, (A7) 



d T pM T ) + \e\d T A(T) • V k p Mk (T)/c 

+ i[H£ t (i) + ^ hoton , r (k, T), p^(T)} 



with p Mk (T) = -iG' i (k,0 - ,T) being the local den- 
sity matrix of electrons with momentum k and 
^photon r(^' ^) coming from the electron-optical- field in- 
teraction given by Eq. ([5]). \e\drA(T)-Vkpnk(T)/c comes 
from the transformation of #tG m (t, T, r) and is just the 
drift term since E(T) = — 9yA(T)/c under the vector 
potential gauge. Moreover, the derivation here is based 
on the assumption that the field A is small and can be 
treated perturbatively. Then, expanding |k| around |k| 
and keeping only the first order term of A(t), one obtains 
Eq. ([6]) except for the scattering term and Hartree-Fock 
term. 



1/ef = l/[e oo + e(q,0)] - l/[e„ + e(q, 0)], (B5) 

with e = 9.7 and e m = 

For the optical phonons, the scattering matrices are 
obtained according to the work by Ishikawa et al^ as 

K^'A = ^W'k/vp^F, (B6) 

with the corresponding parameters of the optical phonons 
given by 



«/i(flq — Ok) 



q- 6 'q)]/2, 



+ ry' e v( 9k +<i _6 'i)]/2, 



AKO 



(B8) 
(B9) 



Appendix B: Expression of the scattering terms 

In this section, we give the scattering matrices used in 
Eqs. d and (HID. C/(e, ;2k -e, ;2k _ q ) = VNlZ l V r (e mk - 
£r; 2 k-q) cxp(— qd) with Zi = 1 being the charge number 
of the impurity, iV, standing for the impurity density and 
d representing the effective distance of the impurity layer 
to the graphene sheet. M£, ? , are scattering matrices 
for electron-phonon scattering with phonons in branch 
A. For the AC phonon, w q Ac = VphQ with v p h being the 



Appendix C: Numerical scheme for drift term 

In this section, we present the numerical scheme for the 
drift term. The truncated 2D momentum space is divided 
into N x M control regions in polar coordinate system,— 
each with equal energy and angle intervals. The grid 
point k n m is localed in the center of the control region 
(n, to) with its radial coordinate k n — (n + 0.5)Afc and 
angular coordinate m = mAQ. Here, Ak is the momen- 
tum interval and AO = 2it /M is the angular interval. 
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The drift term is dealt with the discrete conservation 
principle: 



|e|E(t) • V k p Mk , w | k=k „ m - " T^Tf- 

|e|£7(t) 



kAkAB 



dkdO^ [dk(kp M k tVri i cos6» k ) 



" f.L.nm.rjrj 



kAkA9 



F, 



Afc(sini 



m+0.5 fi^nm^r]' 



sinw m _ .5-f M! „ )m _i^^'JJ 



(CI) 



Here f2 n , OT and 9f2 n)TO are the control region which con- 
tains the grid point k„, m and the corresponding bound- 
ary, respectively, and the electric field is set to be along 
the x direction. In the last step of the above equation, 
the integration of the boundary is replaced by the sum- 
mation over the first order quadrature on the four (or 
three if the control region is the neighbor of k = 0) sides 



of the boundary dQ n ,m an d the fluxes treated up to third 
order~ in the polar coordinate are given by 



%%m,rm' = - fc «'-l cos ^»'/^k„'-i, m ':>)>)7 6 

+ coa9 m (5k n >Pitk. n , im ,T l T l ' + 2fen'+i/Vk„/ +1 . m ,w')/6;(C2) 
ltW '/6. (C3) 



In order to satisfy the requirement of the numerical sta- 
bility, m' in F '™ m , is taken to be m if — |e|E(i) • 



xsmtf m+0 .5 



> and 



1 otherwise. Moreover, n! in 



Funmnri' * s cnosen to be n if — |e|E(i) • xcos# m > and 
n + 1 otherwise. Besides, if n' — 1 > 0, m' is taken to 
be to. On the other hand, if n' — 1 < 0, the term n' — 1 
in the equation is replaced by and to' = to + M/2. It 
is noted that with this numerical scheme, the tempera- 
ture variation due to the noise is within 3 K for the pulse 
investigated in Fig. [2] 
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